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Abstract 

The package fhi98PP allows one to generate norm-conserving pseudopotentials adapted 
to density-functional theory total-energy calculations for a multitude of elements 
throughout the periodic table, including first-row and transition metal elements. 
The package also facilitates a first assessment of the pseudopotentials' transferabil- 
ity, either in semilocal or fully separable form, by means of simple tests carried 
out for the free atom. Various parameterizations of the local-density approxima- 
tion and the generalized gradient approximation for exchange and correlation are 
implemented. 
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PROGRAM SUMMARY 

Title of the program: fhi98PP 
catalogue number: ... 

Program obtainable from: CPC Program 
Library, Queen's University of Belfast, N. 
Ireland (see application form in this issue) 

Licensing provisions: none 

Computers on which the program has been 
developed: IBM/RS 6000, Pentium PC 

Operating system: UNIX 

Graphics software to which output is tai- 
lored: XMGR or XVGR (both are pubhc 
domain packages) 

Programming language: FORTRAN 77 
(non-standard feature is the use of 
END DO); UNIX C-shell scripts are em- 
ployed as command line interfaces 

Memory required to execute with typical 
data: < 1 Mbyte 

No. of bits in a word: 32 

Memory required for test run: < 1 Mbyte 

Time for test run: ~ 1 min 

Number of lines in distributed program: 
7500 (19000 with hnear algebra library) 

PACS codes: 71.15.Hx, 71.15.Mb, 71.15.- 
m, 82.20.Wt 

Keywords: pseudopotential, total energy, 
electronic structure, density functional, lo- 
cal density, generalized gradient 

Nature of the physical problem 

The norm-conserving pseudopotential con- 
cept allows for efficient and accurate ab 
initio electronic structure calculations of 
poly-atomic systems. The key features of 
this approach are: (i) Only the valence 
states need to be calculated. The core 
states are considered as chemically in- 
ert, and removed within the frozen-core 



approximation. This exploits that many 
chemical and physical processes are gov- 
erned by the valence states but only in- 
directly involve the core states, (ii) The 
valence electrons move in a pseudopoten- 
tial which is much smoother than the 
true potential inside the small core re- 
gions around the nuclei, while reproduc- 
ing it outside. This pseudopotential acts 
on smooth pseudo wavefunctions that are 
equivalent to the true valence wavefunc- 
tions, but avoid the radial nodes that 
keep the true valence and core orbitals or- 
thogonal. This enables the use of compu- 
tationally expedient basis sets like plane 
waves, and facilitates the numerical solu- 
tion of the Schrodinger and Poisson equa- 
tions in complicated systems, (iii) The 
norm-conservation constraint ensures that 
outside the core the pseudo wavefunctions 
behave like their all-electron counterparts 
over a wide range of different chemical 
situations. Along with a proper design, 
this makes the pseudopotential approach 
a dependable approximation in describing 
chemical bonds. 

Derived and applied within density- 
functional theory [1-3], norm- conserving 
pseudopotentials [4-6] enable total-energy 
calculations of complex poly-atomic sys- 
tems [7-9] for a multitude of elements 
throughout the periodic table. Questions 
addressed with pseudopotentials provided 
by this code, or its earlier version, range 
from phase transitions [10,11], defects in 
semiconductors [12-14], the structure of 
and diffusion on surfaces of semiconduc- 
tors [15-17], simple metals [18], and tran- 
sition metals [19-21], up to surface reac- 
tions [22,23], including molecules [24,25] of 
first-row species. 

This package is a tool to generate and val- 
idate norm-conserving pseudopotentials, 
usable either in semilocal or in fully 
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separable form, and including relativis- 
tic effects. Exchange and correlation is 
treated in the local-density approximation 
based on Ceperley and Alder's data [26] 
as parametrized, e.g., by Perdew and 
Wang [27], or in the generalized gradient 
approximation, as proposed by Perdew, 
Burke, and Ernzerhof (PBE) [28], Perdew 
and Wang (PW91) [29], Becke and Perdew 
(BP) [30,31], and by Lee, Yang, and Parr 
(BLYP) [32]. 
Method of solution 

The first part of the program (psgen) 
generates pseudopotentials of the 
Hamann [33] or the TrouUier- Martins 
type [34], based on a scalar-relativistic 
all-electron calculation of the free atom. 
A partial core density can be included to 
allow for nonlinear core- valence exchange- 
correlation [35] where needed, e.g., for 
spin-density functional calculations, alkali 
metal compounds, and the cations of II- VI 
compounds like ZnSe. The second part 
(pswatch) serves to assess the transfer- 
ability of the pseudopotentials, examining 
scattering properties, excitation energies, 
and chemical hardness properties of the 
free pseudo atom. Transcribing the pseu- 
dopotentials into the fully separable form 
of Kleinman and Bylander [36], we verify 
the absence of unphysical states by in- 
spection of the bound state spectrum and 
by the analysis of Gonze et al. [37]. The 
convergence of the pseudo wavefunctions 
in momentum space is monitored in order 
to estimate a suitable basis set cutoff in 
plane wave calculations. 
Restrictions on the complexity of the prob- 
lem 

(i) Only some of the GGA's currently in 
use are implemented, others may be read- 
ily added however, (ii) The present pseu- 
dopotentials yield the correct relativistic 
valence levels where spin-orbit splittings 



are averaged over, as it is intended for most 
applications. 
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1 Introduction 



Electronic structure and total-energy calculations using density-functional theory [1-3] 
have grown into a powerful theoretical tool to gain a quantitative understanding of the 
physics and chemistry of complex molecular, liquid, and solid state systems. The norm- 
conserving pseudopotential approach provides an effective and reliable means for per- 
forming such calculations in a wide variety of poly-atomic systems, particularly, though 
not only, together with a plane wave basis and modern minimization algorithms for the 
variational determination of the ground state [7-9,38]. In this approach only the chemi- 
cally active valence electrons are dealt with explicitly. The chemically inert core electrons 
are eliminated within the frozen-core approximation [39], being considered together with 
the nuclei as rigid non-polarizeable ion cores. In turn, all electrostatic and quantum- 
mechanical interactions of the valence electrons with the ion cores (the nuclear Coulomb 
attraction screened by the core electrons, Pauli repulsion and exchange-correlation be- 
tween core and valence electrons) are accounted for by angular momentum dependent 
pseudopotentials. These reproduce the true potential and valence orbitals outside a cho- 
sen core region but remain much weaker inside. The valence electrons are described by 
smooth pseudo orbitals which play the same role as the true orbitals, but avoid the nodal 
structure near the nuclei that keeps core and valence states orthogonal in an all-electron 
framework. The respective Pauli repulsion largely cancels the attractive parts of the true 
potential in the core region, and is built into the therefore rather weak pseudopoten- 
tial [40] . This "pseudoization" of the valence wavefunctions along with the removal of the 
core states eminently facilitates a numerically accurate solution of the Schrodinger and 
Poisson equations, and enables the use of plane waves as an expedient basis set in elec- 
tronic structure calculations. By virtue of the norm-conserving property [4,5] and when 
constructed properly pseudopotentials present a rather marginal approximation, and in- 
deed allow for an adequate description of the valence electrons over the entire chemically 
relevant range of systems, i.e. atoms, molecules, and solids. 

Here we supply a two-piece package to generate norm-conserving pseudopotentials 
adapted to density-functional electronic structure calculations, and usable in semilocal 
or fully separable form. Generic input files for many elements are distributed along with 
the package. 

Part 1 (program psgen) serves to generate the pseudopotentials using the schemes by 
Hamann [33] or by Troullier and Martins [34]. This combination provides efficient pseu- 
dopotentials for "canonical" applications like group IV and III-V semiconductors [9,10,15- 
17,22,23], as well as for systems where first-row, transition or noble metal elements are 
present [11,19-21,24,25], or where "semicore" d-states must be treated as valence states, 
hke in GaN [41] or InN [14]. In such cases the strongly localized 2p and 3, 4, 5d valence 
states are readily handled with the TrouUier-Martins scheme. Both schemes have been 
routinely employed and have an established track record with the electronic structure 
and ab initio molecular dynamics package fhi96md [9]. The pseudopotentials are derived 
within density-functional theory, starting from a scalar-relativistic all-electron calcula- 
tion of the free atom [42] . Accordingly they yield the proper relativistic positions of the 
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valence levels [43] with the spin-orbit coupling being averaged over, in line with the prac- 
tice in most applications. The exchange-correlation interaction may be described within 
the local-density approximation (LDA) [44,45], or the generalized gradient approximation 
(GGA) [46]. For the LDA, the package includes the expressions by Perdew and Wang [27], 
and by Perdew and Zunger [47] , both parameterizing Ceperley and Alder's results for the 
homogeneous electron gas [26]. Earlier prescriptions by Hedin and Lundquist [48], and 
Wigner [49] are supplied as well. For the GGA, the package includes the formulations 
by Perdew, Burke, and Ernzerhof (PBE) [28], by Perdew and Wang (PW91) [29], as 
well as those composed from the exchange functional by Becke [30] (or that of Ref. [29]) 
and the correlation functional either by Perdew (BP) [31], or by Lee, Yang, and Parr 
(BLYP) [32]. Alternative forms can be readily added to the code. Hence the pseudopo- 
tentials, i.e. the electron-ion interaction, can be generated consistently using the same 
exchange- correlation scheme as for calculating the poly-atomic system [50]. A partial core 
density may be included in the unscreening of the pseudopotentials to explicitly account 
for the nonlinearity of core- valence exchange-correlation of the above functionals [35]. 
This has proven necessary, e.g., for alkali metal elements [51], the cations in strongly 
polar compounds like ZnSe [52], or calculations of spin-polarized systems [53]. 
Part 2 (program pswatch) facilitates the assessment of the pseudopotentials' transferabil- 
ity, that is, their proper performance in different chemical environments. This is done 
for the free (pseudo) atom by checking that its scattering properties, excitation energies 
and chemical hardness properties [54,55] faithfully represent those of the correspond- 
ing all-electron atom. For this purpose, we evaluate the logarithmic derivatives of the 
valence wavef unctions, total energies and orbital eigenvalues for different valence config- 
urations. For the commonly used fully separable Kleinman-Bylander form of the pseu- 
dopotentials [36], it is important to rule out "ghost states", i.e. unphysical bound states 
or resonances appearing in the valence spectrum. As an immediate check we evaluate the 
bound state spectrum of the fully separable potentials, and apply the analysis by Gonze 
et al. [37]. To assess the convergence behavior of the pseudopotentials we monitor the 
pseudo atom's kinetic energy in momentum space [56], which also allows a first estimate 
of a suitable basis cutoff for plane wave calculations. 

In Chapter 2 of this Gommunication we briefly survey the conceptual aspects of construct- 
ing and validating pseudopotentials, more practical topics being addressed towards the 
end of a section. Ghapter 3 is meant as a reference guide to our package. 

2 General background and formalism 

A pseudopotential is constructed to replace the atomic all-electron potential such that core 
states are eliminated and the valence electrons are described by nodeless pseudo wave- 
functions. In doing so the principal objectives to consider are (i) the transferability of the 
pseudopotential, its ability to accurately describe the valence electrons in different atomic, 
molecular, and solid-state environments. In self-consistent total energy calculations this 
means that the valence states have the proper energies and lead to a properly normalized 
electron distribution which in turn yields proper electrostatic and exchange-correlation 
potentials, particularly outside the core region, i.e. where chemical bonds build up. (ii) 
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their efficiency, that is to keep the computational workload in applications as low as possi- 
ble, allowing to compute wavcfunctions and electron densities with as few basis functions 
(being "soft" ) and operations as possible. Modern norm-conserving pseudopotentials, im- 
plemented here in the variants proposed by Hamann [33] and by TrouUier- Martins [34], 
allow for a controlled "compromise" with respect to these (conflicting) tasks, through 
appropriate choices and tests at several stages of their actual construction. 

Norm- conserving pseudopotentials are derived from an atomic reference state (Sec. 2.1 
and 2.2), requiring that the pseudo and the all-electron valence eigenstates have the same 
(reference) energies and the same amplitude (and thus density) outside a chosen core cutoff 
radius r'^. Normalizing the pseudo orbitals then implies that they include the same amount 
of charge in the core region as their all-electron counterparts, being "norm-conserving" . 
At the reference energies, the pseudo and all-electron logarithmic derivatives (which play 
the role of boundary conditions on the wavcfunctions inside and outside the core region) 
agree for r > r'^ by construction. Norm-conservation then ensures that the pseudo and 
all-electron logarithmic derivatives agree also around each reference level, to first order in 
the energy, by virtue of the identity 



where ip{e; r) denotes the regular solution of the Schrodinger equation at some energy e 
(not necessarily an eigenfunction) for either the all-electron potential or the pseudopo- 
tential. A norm-conserving pseudopotential thus exhibits the same scattering properties 
(logarithmic derivatives) as the all-electron potential in the neighborhood of the atomic 
eigenvalues, typically it does so over the entire energy range of valence bands or molecular 
orbitals, which is an important prerequisite and measure for the pseudopotential's proper 
performance. 

Regarding transferabiUty the quaUty of a pseudopotential depends on (a) correct scat- 
tering properties, (b) the actual choice of the core cutoff radius r'^, larger values leading 
to softer pseudopotentials (more rapidly convergent in say plane wave calculations) but 
tending to decrease their accuracy, as the pseudo wavcfunctions and the pseudopotential 
begin to deviate from their true counterparts even at radii that are relevant to bonding, (c) 
an adequate account of the nonlinear exchange-correlation interaction between core and 
valence electrons (Sec. 2.3), which is often treated hnearly through the pseudopotential 
and thus approximated, (d) if applied, the quality of the Kleinman-Bylander form of the 
pseudopotential, where accidental spurious states need to be avoided (Sec. 2.6), (e) the 
validity of the frozen-core approximation, i.e. the inclusion of all relevant states as valence 
states which may include upper semicore states in addition to the outermost atomic shell, 
dependent on the actual chemical context, say, if these hybridize with valence states of 
neighboring atoms (Sec. 2.3), and (f) the treatment of higher angular momentum compo- 
nents which are just implicitly considered in the pseudopotential's construction (Sec. 2.5). 
As outlined in Sec. 2.4 the points (a) - (e) can be controlled with respect to transferability 
by verifying that the pseudo atom reproduces the all-electron atom's scattering properties 
and excited states. The latter essentially probe the atom's response or hardness properties 
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due to a change in the occupation of the valence orbitals [54]. Such tests for the free atom 
are less conclusive for (c) and (f) which, though mostly uncritical, should be kept in mind 
as potential sources of pseudopotential errors in poly-atomic systems. 

Usually norm-conserving pseudopotentials are set up and derived in terms of an angular 
momentum dependent, "semilocal" operator, 

'max m=l r/ _ n 

(r|yP^|r')=y^-(r)5(r-rO + E yL{^r)SVr{r)^-^YU^r,) , (2) 

1=0 m=-l ' 

written in terms of a local pseudopotential V^°'^{r), and /-dependent components SVi^^{r) ~ 
VJ^^(r) — V^°^{r) which are confined to the core region and eventually vanish beyond some 
/max (see Sec. 2.5). For reasons of computational efficiency the semilocal form is often 
transcribed into the fully nonlocal Kleinman-Bylander form [36] (Sec. 2.6). Relativistic 
effects on the electrons, important for heavier elements, basically originate just in the core 
region. Whence they may be incorporated in the pseudopotentials [57,58] which we do 
here within a scalar- relativistic approximation (Sec. 2.1). Accordingly, such "relativistic" 
pseudopotentials are formally employed with a non-relativistic Schrodinger equation. 
The next subsections discuss the construction of this package's pseudopotentials, which 
in short proceeds as follows, 

• Density-functional calculation of the all-electron atom in a reference state (ground 
state), within a scalar- relativistic framework and a chosen approximation for exchange- 
correlation (Sec. 2.1). 

• Construction of the pseudo valence orbitals and (screened) pseudopotential components, 

observing the norm-conserving constraints (Sec. 2.2). 

• Removal of the electrostatic and exchange-correlation components due to the valence 
electrons from the screened pseudopotential, optionally taking nonlinear core-valence 
effects into account. This "unscreening" yields the ionic pseudopotential representing 
the electron- ion interaction in poly-atomic systems (Sec. 2.3). 

• Assessment of the pseudopotential's transferability (Sec. 2.4), in case of a transforma- 
tion to the fully separable Kleinman-Bylander form the exclusion of unphysical "ghost 
states" in the valence spectrum (Sec. 2.6). 

Hartree atomic units are used throughout {e"^ = h = m = 1), lengths being given in units 
of the Bohr radius 1 bohr = 0.52918 A, and energies in units of 1 hartree = 27.2116 eV. 

2.1 Atomic all-electron calculation 

The initial step in constructing our pseudopotentials is an all-electron calculation of the 
free atom in a reference configuration, usually its neutral ground-state. Using density- 
functional theory to handle the many-electron interactions, the total energy of N electrons 
in an external potential V^^^{r) (for an atom the nuclear —Z/r potential) is written as a 
functional of the electron density p, 

E'^'ip] = T[p] + £;^C[p] + E'^ip] + J ~p{v)dh , (3) 
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where T[p] is the non-interacting kinetic energy, E'^'^[p] the exchange and correlation en- 
ergy, and = I / ^ |r^-r'"| ^ '^'^^ ^'^^ electrostatic or Hartree energy from the electron- 
electron repulsion. The ground-state corresponds to the minimum of E^°^[p], and is found 
from the self-consistent solution of the Kohn-Sham equations 



--V^ + y^"[p;r]-e, 



^,(r) = 



(4) 



Sp{v) 



(5) 



P(r) = E/^l^^(r)r , jp{v)dh^N , 



(6) 



fi^l if ej < ejv , < /i < 1 if = , and /i = if > e^r . (7) 

Specializing to a spherically symmetric electron density and effective potential the Kohn- 
Sham orbitals separate, ilJiir) = [M„.;.(ei; r)/r]y;.mi(^^r) and (4) reduces to a radial 
Schrodinger equation. Relativistic effects, which in principle require a four-current and 
Dirac-spinor formulation of Eqs. (3) - (7), are treated by using a scalar- relativistic ki- 
netic energy operator [42]. This takes into account the kinematic relativistic terms (the 
mass-velocity and Darwin terms) needed to properly describe the core states and the rel- 
ativistic shifts of the valence levels, in particular for the heavier elements. The spin-orbit 
coupling terms are averaged over, as is commonly done in most applications (see also 
Ref. [58]). The radial wavefunctions Uni are then the self-consistent eigenfunctions of the 
scalar-relativistic Schrodinger equation 



1 dV{r) d 1 /(/ + 1)\ , 



2M(r) V dr'^ 2M{r)c^ dr drr r 



Uni{^i;r)^0, (8) 



with the relativistic electron mass M(r) = 1 + ^^^1^, where 1/c = 1/137.036 is the 
finestructure constant, the effective all-electron potential 

V{r) = V''''[p;r] + V''\p;r]-^ , (9) 
the exchange-correlation potential l^^*"[p;r] = ^^p(^r)'^ , and the Hartree potential 

V'^lp; r] = = 47r- p{ry^dr' + Att H p{r')r'dr' . (10) 

Opyr) r Jo Jr 

Due to the spherical symmetry, states with the same quantum numbers Ui = n and 
li — I , but different are energetically degenerate, = emkmi = ^ni, and therefore 
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equally occupied by < /j < fni/{'2l + 1) electrons, where the occupancies /„; stand for 
the total number of electrons in the nZ-shell. This leads to the electron density 

fni = 2(2/ + 1) for e„/ < cn , 
0<fni< 2(21 + 1) for e„; = ejv, (H) 
/„; = for €nl > ejv , 

with the total number of electrons being J2ni fni = N . The non-relativistic Kohn-Sham 
equations are recovered by setting the finestructure constant 1/c in Eq. (8) to zero. The 
above scheme determines the atomic ground state. Excited atomic states are calculated 
with appropriately specified occupancies (as usual, we treat the occupancies always as 
parameters, a database of ground state configurations {fni} being supplied with this 
package). 

Regarding exchange-correlation, the program allows to employ commonly used param- 
eterizations of the local-density approximation (LDA) and of the generalized gradient 
approximation (GGA). In the LDA the exchange-correlation energy is 



n— 1 



n=l 1=0 



, with 



^XC-LDA[^] ^ I ^XC-LDA(^(^)y3^ ^ (12) 

where e^^'~^^^(p) = e^~^^^{p) + e'~"~^^^{p) represents the exchange-correlation energy 
per electron of a uniform electron gas. Its exchange part reads s^^^^^{p) = — J;:(37r^p)^/^. 
For the correlation part e^~~^^^{p) representations given by Pcrdew and Wang [29] and, 
previously, by Perdew and Zunger [47] may be used. Both arc parameterizations of Ceper- 
ley and Alder's exact results for the uniform electron gas [26] . To facilitate comparisons 
with results in the existing literature, we also supply the earlier prescriptions by Hedin 
and Lundquist [48], and by Wigner [49]. For an overview of the LDA wc refer to Rcf. [45]. 
These (nonrelativistic) LDA exchange-correlation functional may be supplemented with 
relativistic corrections to the exchange part, as given by MacDonald and Vosko [59,60]. 
The LDA exchange-correlation potential is given by 



V 



XC-LDA 



[p; r] 



1 + p(r) 



dp 



e 



XC-LDA 



iPir)) 



(13) 



In a GGA the exchange-correlation functional depends on the density and its gradient, 

^XC-GGA[^] _ I ^(^^^XC-GGA(^(^)^ |Vp(r)|) . (14) 

This code includes the GGA functionals by Perdew and Wang (PW91) [29], by Perdew, 
Burke, and Ernzerhof (PBE) [28], and Becke's formula for the exchange part [30] combined 
with Perdew's 1986 formula for correlation [31] (BP). Substituting for the latter the 
formula of Lee, Yang, and Parr (LYP) [32] provides the so-called BLYP GGA. Relativistic 
corrections to the GGA as proposed by Engel et al. [61] are not included here. The GGA 
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exchange-correlation potential (in cartesian coordinates) is given by 



p(r)£^^-^^^(p(r),|Vp(r)|) , (15) 

and depends on the first and second radial derivatives of the density. It is an open issue 
whether to prefer one of these GGAs over another. Dependent on the application, they 
may yield somewhat differing results. However, the actual choice among these GGAs is 
usually less important than the differences between the LDA and the GGAs themselves. 
A discussion of the various GGAs can be found in Ref . [46] . 

Whether to use either the LDA or the GGA in constructing pseudopotentials is determined 
by the exchange-correlation scheme one wants to use for the poly-atomic system. If using 
the GGA there, the pseudopotential should be generated within the (same) GGA rather 
than the LDA, in order to describe the exchange-correlation contribution to the core- 
valence interactions within the GGA too. This "consistent" use of the GGA may be 
significant to resolve any differences between LDA and GGA properly, i.e. to the same 
extent as with all-electron methods [50] (see also Sec. 2.3). 



V 



XC-GGA 



d ^ d 
dp fr{ dViP 



2.2 Screened norm- conserving pseudopotentials 



Having obtained the all-electron potential and valence states, we use either the scheme 
by Hamann [33], or by Troullier and Martins [34] to construct the intermediate, so called 
screened pseudopotential. The latter acts as the effective potential on the (pseudo) valence 
states ipf^ir) = [uf^{^f^',r)/r] Yi^^flj.), where the radial part for the angular momentum / 
is the lowest eigenfunction of the nonrelativistic Schrodinger equation. 



1 l{l + 1) 



+ Vr-{T) - e\ 



(16) 



Initially, a radial pseudo wavefunction uY{r) is derived from the all-electron valence level 
with angular momentum / (the reference level [62]) such that 

(1) The pseudo wavefunction and the all-electron wavefunction correspond to the same 
eigenvalue 



(17) 



and their logarithmic derivatives (and thus the respective potentials) agree beyond 
a chosen core cutoff radius rf , 



^ In uf{ef] r) ^ ^ In Uni{enh r) for r > r\ 



(18) 



(2) The radial pseudo wavefunction has the same amplitude as the all-electron wave- 
function beyond the core cutoff radius. 



^r(er; Uni{eni] r) for r > r\ 



(19) 
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and is normalized, \uf^{ef^; r)\^dr — /q°° \uni{eni', r)pdr = 1. This implies the norm- 
conservation constraint 

r \uf{er;r)\'dr^ r Ki{eni;r)\'dr , for r' > r[ . (20) 
Jo Jo 

(3) The pseudo wavefunction contains no radial nodes. In order to obtain a continuous 
pseudopotential which is regular at the origin, it is twice differentiable and satisfies 

The pseudopotential components then correspond to an inversion of the Schrodinger equa- 
tion (16) for the respective pseudo wavef unctions, 

v;-M^er-^ + ^^»rM . (21) 

and become identical to the all-electron potential beyond r > rf. 

The ansatz for the spatial function uf^{r) at the energy ef^ requires three free parameters 
to meet the constraints (i) and (ii). The Hamann scheme is of this "minimal" type, 
while the Troullier-Martins scheme introduces additional constraints, namely that the 
pseudopotential's curvature vanish at the origin, ^Vi^''"'\r)\r=o = 0, and that all first 
four derivatives of the pseudo and the all-electron wavefunction agree at the core cutoff 
radius. Thereby it achieves softer pseudopotentials for the 2p, and the 3,4, Sc? valence 
states of the first row, and of the transition metal elements, respectively. For other elements 
both schemes perform rather alike. We refer to the original references [33] and [34] for the 
details about the construction procedure. Here we point out that both schemes technically 
differ somewhat regarding the radii where the pseudo wavefunctions and pseudopotentials 
match their all-electron counterparts. The Hamann scheme accomplishes the matching 
exponentially beyond the core cutoff radius r'^~^, while in the Troullier-Martins scheme 
this matching is exact at and beyond r[~^^^. Although in the Hamann scheme the core 
cutoff radii are nominally smaller, r^^^/r^^™ 0.25 . . .0.75, the pseudo wavefunctions 
and pseudopotentials converge to the all-electron wavefunctions and potential within a 
similar range as in the Troullier-Martins scheme. 

Default values for the core cutoff radii are provided by the program psgen. In the Hamann 
case it determines the position r^^^ of the outermost maximum of the all-electron wave- 
function and sets = 0.6 rf'^ if there are core states present with the same angular 
momentum, and to rp^ = 0.4 rf^^ otherwise [33] . In the Troullier-Martins case the core 
cutoff radii are set to values which we have tabulated for all elements up to Ag (except 
the lanthanidcs). Except for some finetuning that might be necessary when the pseudopo- 
tentials are used in the Kleinman-Bylander form (see Sec. 2.6) these default values should 
yield a reasonable compromise between pseudopotential transferability and efficiency. 
In general, increasing r[ yields softer pseudopotentials, which converge more rapidly e.g. 
in a plane wave basis. These will become less transferable however as the pseudo wave- 
functions become less accurate at radii relevant to bonding. When decreasing rf one has 
to keep in mind that rf has to be larger than the radius of the outermost radial node 
in the respective all-electron orbital. Taking rf too close to such a node results in poor 
pseudopotentials, as it is essentially at variance with requiring both a nodeless and at 
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the same time norm-conserving pscudo wavcfunction. A rugged or multiple- well structure 
seen in a plot of the screened pseudopotential and poor scattering properties (see Sec. 2.4) 
can indicate such a breakdown. 

Treating components without bound reference states 

So far we have considered pseudopotentials for bound reference states. On the other 
hand it is quite common in LDA and GGA that unoccupied states in the neutral atom 
are not or only weakly bound, for instance the 3(i-wavefunctions for Al, Si, etc. . To 
derive a /-component of the pseudopotential in such cases we follow the generalized norm- 
conserving procedure proposed by Hamann [33]: At a chosen energy ei the Schrodinger 
equations for the all-electron potential and the pseudopotential are integrated from the 
origin up to a matching radius r™'^*'^'^. There the logarithmic derivatives and the ampli- 
tudes of the pseudo and all-electron partial waves are matched, and the "normalization" 

^. match ^match 

Jo' \uY {^h r)^^ dr = Jq' \uni{ef,r)\'^dr := 1 is enforced. This procedure is analogous 
to (17) - (20) for bound states. The corresponding pseudopotential component results 
from (21), and merges into the all-electron potential for r > rf^^^^^. The value of €i should 
be chosen in the energy range where the valence states are expected to form bands or 
molecular orbitals, as a default, the program psgen employs the highest eigenvalue of 
the occupied states, e.g., for Al to e;=2 = ^3p- Hamann's procedure for unbound states 
is not tied to a particular type of pseudopotential, we use it in constructing Troulher- 
Martins as well as Hamann pseudopotential functions. Figure 1 displays the Hamann-type 
pseudo wavefunctions and pseudopotential components for aluminum. Alternatively, one 
can derive the /-components without bound states in the neutral atom from a (separate) 
calculation of an ionized or excited atom which supports bound states for these /'s [6]. 
The above prescriptions should both yield rather similar pseudopotentials, consistent with 
the supposition that pseudopotentials are transferable among neutral and excited atomic 
states. 

2.3 Ionic pseudopotentials, unscreening, and nonlinear core-valence exchange- correlation 

The final ionic pseudopotentials are determined by subtracting from the screened pseu- 
dopotentials the electrostatic and the exchange-correlation screening contributions due to 
the valence electrons, 

Vr{r) = Vr'\r) - r] - V^'^l^- r] with pl\r) = ^Y.fi 

^'^ 1=0 

(22) 

The valence electron density is evaluated from the atomic pseudo wavefunctions, with the 
same occupancies fi as for the all-electron valence states of Sec. 2.1. Figure 1 shows the 
ionic pseudopotentials for aluminum as an example. With the pseudopotentials as defined 
by Eq. (22) the electronic total energy of an arbitrary system reads 

E'°'^YlM^i\'^ + V'-'\A) + E''[pn+E''''[pP^] with pP^(r) = X:/#i(r)P . (23) 

i i 
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Pseudo vs All-Electron Wavefunctions Al Ionic Pseudopotentials A! 

Fig. 1. Graphical output from psgen for aluminum. The left panel shows the valence wave- 

j 'unctions for the i pseudo atomi (sohdi lines) and tho aU-ologtrpn latQWii (dashed). i The i right panel 



o 

S-l 

cS 

!-l 
4J 
• 1— I 

!-l 

cS 



shows the respective Hamann-type ionic pseudopotential, using Hamann's procedure to treat 

give the core cutoff radii r[ (in bohr^. ~ ^ 



the d-compon^t (lee Wc?. 2.2) and the LDA. The 

3s 1 



1.0 



2p 
3p 
3d 
3d 



1.547 



1.547 



0.0 



212 




Here, E-^'^ refers to the exchange-correlation interaction between the valence electrons 
themselves. That between the valence and the core electrons is included in the pseudopo- 
tential energy, as a term which depends linearly on the valence density p^^. Although E^'^ 
is a nonlinear functional of the total electron density p, the above "linearization" of its 
core-valence contribution is a usual and mostly adequate approximation for calculations 
within both LDA and GGA [50]. However, an explicit account of the core- valence nonlin- 
earity of E^^ is sometimes required, for instance in studies involving alkali metal atoms or 
within spin-density functional theory. This is accomplished by restoring the dependence 
of E^'^ (and V^'^) on the total electron density, i.e. adding the (atomic) core density 
to the valence density in the argument of E^'^. In turn the ionic pseudopotentials (22) 
are redefined. Rather than the full core density, it suffices to add a so called partial core 
density Po""^^, as suggested by Louie et al. [35]. It reproduces the full core density p^'^^ 
(47) outside a chosen cutoff radius r°'^ but is a smoother function inside, which enables 
its later use together with plane waves. We use a polynomial 



prir) = 



for r > r 



nlc 



Co + Ei=3 CiT' for r < r''^'' , with p^'^r) < p^'^r) , 



(24) 
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where we take the coefficients q such that pQ°^'^{r) has zero slope and curvature at the 
origin, decays monotonically, and joins the full core density continuously up to the third 
derivative. In this way smooth exchange-correlation potentials result for the LDA as well 
as for the GGAs. The resulting nonhnear core-valence exchange-correlation scheme uses 
the redefined ionic pseudopotentials 

Vrir) ^ Vr'^'ir) - r] - V'Vo + Pr ; r] (25) 

instead of Eq. (22). The according total energy expression is 

E""' ^ E M'^ilT + V'^'IA) + E"[p^'] + + pri , (26) 

i 

where, in a poly-atomic system, p^^"^ stands for the superimposed partial core densities 
of the constituent atoms. 

Whether nonlinear exchange-correlation plays a significant role can be readily decided 
by a comparison of two calculations with and without a partial core density. As a rule, 
it becomes more important to the left in the periodic table, i.e. with fewer electrons in 
the valence shell (alkali metals), and with increasing atomic number, i.e. the farther the 
upper core orbitals extend into the tail of the valence density (e.g. Zn, Cd). Of course 
if the uppermost semicore states hybridize with "true" valence states they have to be 
treated as valence states in the first place, like the 3, 4, 5(i-states in all transition and 
noble metals, but also the Sd-states in Ca. This also applies for the 3d- and id-states of 
Ga and In in GaN and InN, where these interact markedly with the N atoms' 2s-states in 
GaN and InN, while treating them in the core is indeed adequate in, say, GaAs or In As. 

There remains the question what value to pick for the core cutoff radius r^^^. It is reason- 
able to set r^^'^ to about the radius where the full core density drops below the valence 
electron density. The utility psgen determines this "equi-density radius" and also outputs 
the core and valence densities. Choosing r"''^ smaller makes it more demanding to represent 
the partial core, say, in a plane wave basis, without further enhancing the pseudopoten- 
tial's performance. On the other hand, too large a value for r^^^ might be insufficient to 
capture the nonlinearities in the desired way. Of course, a systematic convergence test can 
be carried out by inspecting the variation with respect to r"^'*^ of some simple system's 
properties, e.g., the crystal equilibrium volume. 

The GGA, unlike the LDA, occasionally gives rise to short-ranged oscillations in the ionic 
pseudopotentials, particularly when Eq. (22) is employed together with the PW91 GGA, 
but less so when a partial core density is included. These features may be traced to the 
behavior of V^^~^^^[pq^; r] subtracted out in the unscreening, and are certainly artifacts 
of the GGA. From the viewpoint of plane wave calculations they correspond to an energy 
regime where the kinetic energy dominates all other contributions and are, therefore, 
physically unimportant. A "controlled" smoothing is then implied already by the plane 
wave basis cutoff which, in our experience, can be chosen very similar to that for the LDA. 
Accordingly we do not apply any smoothing of the GGA pseudopotentials in real space. 
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2-4 Transferability and convergence behavior 



Transferability considerations - By construction, a pseudopotential reproduces the va- 
lence states of the free atom in the reference configuration. In apphcations, it needs to 
perform correctly in a wide range of different environments, say free molecules and bulk 
crystals, "predicting" the same valence electronic structure and total energy differences 
that an all-electron method would deliver. This transferability depends critically on (a) 
the choice of the core cutoff radii (see Sec. 2.2), (b) the linearization of core- valence 
exchange-correlation (Sec. 2.3), (c) the frozen-core approximation underlying the pseu- 
dopotential's construction (Sec. 2.3), and (d) the transformation of the semilocal into the 
fully separable form of the pseudopotential (discussed in Sec. 2.6). 

Clearly, the transferability of the pseudopotentials should be carefully verified before em- 
barking on extensive computations. Below we discuss several simple tests to be performed 
for a free atom. These tests assume spherical symmetry, and do not mimick non-spherical 
environements. It is therefore good practice to also cross check - for some properties of 
simple test systems, e.g., bond lengths in molecules or a bulk crystals - the results from 
one's pseudopotential calculation with all-clcctron results, which often are available from 
the literature. The compared data have to refer to the same approximation for exchange- 
correlation, of course. Note that a cross check against experimental data alone might be 
less conclusive, since the LDA or GGA themselves lead to systematic deviations between 
theoretical and experimental figures which need to be distinguished from errors rooted in 
the pseudopotentials themselves. 

Test of scattering properties - As a first test, we establish that the logarithmic derivatives 
of the radial wavef unctions. 



A(6,r^^'^^)= ^]nui{e;r 
ar 



(27) 



«diag 



agree for the pseudo and the all-electron atom, as a function of the energy at some diag- 
nostic radius r^'^^ outside the core region, say half of a typical interatomic distance. The 
norm-conservation feature ensures that the logarithmic derivatives for the pseudopoten- 
tial (semilocal or fully separable) are correct to first order around the reference energies 



, by Eqs. (1) and (20), 



d d 
de dr 



liiUni{e;r) 



d d pg, , 
d-ed'r^''''^ ^^'^^ 



, with r > r[ 



(28) 



More generally the pseudopotential and all-electron logarithmic derivatives should closely 
agree over the whole range of energies where the valence states are expected to form 
Bloch bands or molecular orbitals, say ±1 hartree around the atomic valence eigenvalues. 
Figure 2 shows, as a typical case, the logarithmic derivatives for an aluminum pseudopo- 
tential, evaluated with this package's utility pswatch. Note that the logarithmic derivatives 
are evaluated for the screened potentials of the atomic reference configuration, and will 
look different for an excited configuration. 

Test of excitation energies - Next we verify that the pseudopotential reproduces the all- 
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(29) 



where {/,^'"} denote the orbital occupancies in the excited and in the ground state config- 
uration, respectively. The various expressions for E^°^~^ are collected in the Appendix. 
This test seeks to emulate the cost in total energy associated with the orbital hybridization 
in different environments ("promotion energies"). The errors due to the use of pseudopo- 
tentials, AE^^ = E]^^ - E^^ should be compared to AEf^ = Ef^ - E^^, the errors 
that result from an all-electron calculation within the frozen-core approximation. There 
one recalculates just the valence states (with the full nodal structure), but keeps the core 
density fixed as in the ground state configuration {/f }. Due to the variational principle, 
frozen-core calculations yield higher total energy values than unrestricted all-electron cal- 
culations where the core is allowed to relax, £'*°*~^^({/f }) > _E'*°*~'^^({/j^}), the equality 
holding only for {/f }. In this test a transferable pseudopotential should perform with the 
same accuracy as the frozen-core all-electron calculation, |Ai?^^| ~ |Ai?^*-^|, where the 
errors are typically of opposite sign for linearized core-valence exchange-correlation, and 
essentially the same if nonlinear core- valence exchange- correlation is taken into account. 
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Fig. 3. Convergence of the total energy in pseudopotential plane wave calculations of bulk GaAs, 
diamond, and fcc-Cu. Solid lines show the absolute total energy error (per electron) expected 
from the atomic pseudo wavefunctions using the estimate Eq. (32). Open circles show the actual 
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plane wave cutoff energy (Ry) 

Typically, say for the first ionization potential, the errors do not exceed a few ten meV. 
Figure 2 shows the excitation energy error for an aluminum pseudopotential. 

Test of level changes (chemical hardness) - As a last point we verify that the orbital 
eigenvalues for the pseudo atom closely track those of the all-electron atom when the or- 
bital configuration changes. This test provides a discrete sampling of the atomic hardness 
matrix, H^j :— ^^j{{fk}) — df.gf. ^^^iifk}) which ultimately governs the changes of the 
one-particle levels [54,55] and, by Janak's theorem [63], of the total energy as functions 
of the orbital occupancies. Like for the excitation energies, the eigenvalues from the pseu- 
dopotential calculation are expected to agree with the all-electron results to within the 
error found for the frozen-core approximation, as seen in Fig. 2. 



Convergence behavior 

Checking the convergence behavior of the pseudopotentials with respect to the basis size 
is a basic task in calculating realistic systems, where the goal is to achieve convergence 
in total energy differences, e.g. binding energies, rather than the total energy itself. To 
provide some guidance in picking the smoothest out of a set of some similarly transferable 
pseudopotentials, and in choosing an initial estimate for a plane wave basis, the program 
pswatch evaluates the momentum space convergence of the atomic pseudo wavefunctions' 
kinetic energy. This provides a measure for the (absolute) total energy convergence in 
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applications [56]. Using the Fourier decomposition of the pseudo wavef unctions, 

= r kTji{kT)uT{T) dr , (30) 

V TT JO 

consider the error of the momentum space value of the kinetic energy as a function of the 
plane wave cutoff energy E'P" corresponding to the cutoff momentum fcP™ = y/2EP^, 



A£;r(0= / ^\uT{k)\'dk- ufir) 
Jo 2 Jo 



/(/ + !) 



uf{r)dr . (31) 



The program pswatch computes the cutoff energies to achieve AE']^™ < 100, 10, 1 meV. 

For real systems a converged calculation typically corresponds to a plane wave cutoff that 
meets the criterion Ai?]^'°(A;P™) < 100 nieV for the pseudo atom [64]. A corresponding 
estimate of the (absolute) convergence of the total energy in poly-atomic systems is then 
given by the average of the errors AE'J^'" over all states and atoms. 



atom 1 



atom 2 



with the orbital weights Wi. Figure 3 illustrates this estimate. 

2.5 Usage of pseudopotentials: Separation of long-ranged local and short-ranged nonlocal 
components 

Operating with the ionic pseudopotential on the wavefunctions in a poly-atomic system 
requires their projection onto an angular momentum expansion, the second term of Eq. 
(2), which causes the predominant computational effort in calculating larger systems. As 
all components Vi^^{r) at large r reduce to the ionic Coulomb potential, —Z^°^/r, getting 
independent of L It is thus expedient to express the pseudopotential as a multiplicative 
local potential plus only a few Z-dependent, short-ranged "corrections" for I < Imax, 

(r|V^P^|r') = {r\V^°' + SV^y) (33) 

- VZrnr - r') + E E Y:jn^)5Vr\r)^-^YU^^,), (34) 

l=Om=-l ' 

where the local potential is taken as one of the semilocal pseudopotential components, 
V^°'^(r) = Vi^i^ (r), and 6Vf^{r) = V^^^(r) — (r) vanishes beyond rf. Typically one 
truncates the expansion at /max = 2, whereby the core effects on the valence states with the 
same symmetries as the core states are accounted for. The number of projections is reduced 
most if one picks the Zioc = /max component as local potential V^°^{r) = Vi^^^{r) = Vi^^^^ir)- 
Most preferably koc = ^max = 2, the common practice for sp-type materials like silicon. 
We note that the local component of the pseudopotential in principle has to reproduce 
the scattering properties in the higher angular momentum channels with / > /max- Further 
aspects on the choice of the local potential come into play when one uses the fully separable 
form of the pseudopotential as discussed in the next Section. 



18 



2. 6 Fully separable pseudopotentials and ghost states 

Actual electronic structure codes mostly use the ionic pseudopotential in the fully sepa- 
rable form of Kleinman and Bylander (KB) [36] 

{r\VPy) = (r|V'°^ + 5V^^^|r') (35) 
- CW<^(r-0 + E E(r|xL")^r(xL"|r') , (36) 

1=0 m=-l 

where the short-ranged second term is a fully nonlocal rather than a semilocal operator in 
r-space. Starting from a semilocal pseudopotential, the corresponding KB-pseudopotential 
is constructed such that it yields exactly the same pseudo wavefunctions and energies. This 
is accomplished with the projector functions 

(r|xL") = lxi{r)YU^r) = l^^^^^ril ^U^r) , (37) 



and the KB-energies 



|L,PS\T/sl|| 



with |K'5V^''|| = uf{r)5Vi-\rfuf{r)dr [65]. The KB-energy measures the strength 
of a nonlocal component relative to the local part of the pseudopotential, where the 
sign is determined by the KB-cosine —1 < {ipim\')(dm) < 1- Note that X;(r) ^ as 
5V{^^{r) outside the core region. In poly-atomic systems the KB-form drastically 
reduces the effort of operating with the pseudopotentials on the wavefunctions, which 
causes the dominant computational workload in particular for larger systems: for an N- 
dimensional basis (IGi)}, the semilocal form requires the evaluation and storage of ~ 
(iV^ -|- N)/2 matrix elements {Gi\5Vf-\Gj) ^ with the KB-form these factorize, requiring 
just ~ projections (Gilx™) and simple multiphcations (G'i|x™)£'™(x™|G'j). 
In transforming a semilocal to the corresponding KB-pseudopotential one needs to make 
sure that the KB-form does not lead to unphysical "ghost" states at energies below or 
near those of the physical valence states as these would undermine its transferability. Such 
spurious states can occur for specific (unfavorable) choices of the underlying semilocal and 
local pseudopotentials. They are an artifact of the KB-form's nonlocality by which the 
nodeless (reference) pseudo wavefunctions need not be the lowest possible eigenstate, 
unlike for the semilocal form [37]. Formally the KB-form can be generalized to a series 
expansion of a fully nonlocal pseudopotential which avoids ghost states by projecting onto 
additional reference states in Eq. (36) [66], yet at the cost of somewhat increasing the 
computational labor. In practice transferable ghost-frcc KB-pseudopotentials are readily 
obtained by a proper choice of the local component (/loc) and the core cutoff radii (rf) in 
the basic semilocal pseudopotentials. In the next paragraphs we discuss how to identify 
and remove ghost states using the utility pswatch. 
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First of all a ghost state is indicated by a marked deviation of the logarithmic deriva- 
tives of the KB-pscudopotential from those of the respective semilocal pseudopotential 
or all-electron potential. Note that such features can be delicate to detect if the energy 
sampling is too coarse to resolve them or if they occur below the inspected energy interval. 
The program pswatch provides the logarithmic derivatives on an adaptive energy mesh, 
integrating the nonlocal radial Schrodinger equation for the KB-pseudopotential, 

r) + Xi{r)Ei / Xi{r')uf{e, r') dr' ^ , (39) 
Jo 

with El = Ef^^, from the origin up to a diagnostic radius, as described in Ref. [37]. 
Secondly, ghost states are revealed in a comparison of the atomic bound state spectra 
for the semilocal and the KB-pseudopotentials [67], i.e., as extra levels not present in the 
true valence spectrum. The program pswatch solves the respective eigenvalue problem, 
Eq. (39), using Laguerre polynomials as a numerically convenient basis [68]. 
Thirdly, ghost states below the valence states are identified by a rigorous criterion by 
Gonze et al. [37]. The program pswatch evaluates this criterion, here we summarize its 
key points. For a given angular momentum, consider the spectrum {ti{Ei)} of the non- 
local Hamiltonian associated with the KB- form, Eq. (39), as a function of the strength 
parameter Ei. Comparing to the spectrum resulting from the local part of the (screened) 
pseudopotential alone, {ei{Ei = 0)}, the eigenvalues always satisfy the sequences [37], 

io{0)<io{Ei)<ii{0)<hiEi)<... if Ei>0 , (40) 

ioiEi) < eo(0) < iiiEi) < ei(0) < ... if Ei < . (41) 

Now the actual KB-pseudopotential gives Ei = Ef^. The reference eigenvalue ef^ which 
it reproduces by construction corresponds to one of the levels {eo,i,...(-E'i = Ef^^)}. If it 
corresponds to the ground state level, i.e. ef^ = eo{Ei = Ef^^), obviously no lower (ghost) 
state exists. Whether this holds true is decided from (40) or (41) with the help of the 
eigenvalues for the local potential: There is no ghost level beneath the reference level ef^ 
for 

if the reference level lies below the first excited level of the local (42) 
potential, eo(0) < < ei(0) in Eq. (40), 

if the reference level lies below the ground state level of the (43) 
local potential, ef < eo(0) Eq. (41). 

The program pswatch reports ghost states based on this criterion. 

As a practical rule ghost-free KB-pseudopotentials may be obtained with the d-component 
as the local potential except for "two-shell" situations like the transition metal elements 
where typically the s-component is used as the local potential. To eliminate a ghost 
state occurring for some / one may (i) change Zioc, i.e. use a different component of the 
semilocal pseudopotential as the local potential, (ii) adjust the core cutoff radii rf of the 
offending component Vi^^{r) or the local potential Vi^^^{r), respectively. In doing so the 
basic transferability ought to be maintained (see Sec. 2.4). Also, it is expedient to take 
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Table 1 

Analysis of ghost states for the Kleinman-Bylander pseudopotential of Cu employing Eq. (42) 
and (43). The eigenvalus for the local potential are underlined. The semilocal pseudopotentials 
were constructed for the (3(i^° 4s^ 4p°)-configuration, using the Troullier-Martins scheme with 
cutoff radii r^g = rf^2 — 2.1 bohr and r^p ~ 2.3 bohr. Exchange-correlation was treated within 
the LDA. 

l\oc I Ef^^ (eV) bound state energies (eV) ghost? 

2 (3d) 0(4s) 319.67 -219.06 -72.98 -15.95 -4.86 yes 

6o(0) eo(£;f^) ei(0) h{E^^) ^ ef. 



1 (4p) 223.44 -108.75 -0.77 -0.46 no 



eo(0) eo(^r) ^ 41 e~i(0) 



0(4s) 2 (3d) -273.80 -5.33 no 



1 (4p) 51.64 -0.92 -0.77 no 



eo(0) eo(^r) ^ 



^loc = ^max in order to keep the number of the (computationally intensive) projections 
implied by Eq. (36) small. Modifying the above entities seeks to adjust Ef^^^ i.e. the 
strength of the KB-pseudopotcntial relative to the local potential such that either (42) or 
(43) is satisfied. Note that a repulsive enough local potential (achieved e.g. by increasing 
rfj^^ or taking li^c = 0) rules out any ghost state: Such a potential eventually leads to 
a negative KB-energy (38) as 5Vi^\r) = V^\r) - V^°%r) < 0. Also, eo(0) since 
such a local potential gives at most a weakly bound ground state. As Ef-^ < one has 
eP < eo(0) ^ by virtue of (41), ruling out any ghost level. 

Table 1 illustrates the foregoing for a copper pseudopotential (see also Ref. [37] and the 
examples in this package's distribution): the 3, 4, 5(i-component taken as local potential 
^yioc _ acts as a strongly attractive potential on the s-states. Accordingly the 

eigenvalue sequence starts at quite low an energy. In the case of copper the reference 
energy of the 4s-pseudopotential lies already above the first excited level for the local 
potential, ei(0) < e^s- By Eqs. (41) or (40) the Cu 4s-lcvcl corresponds to an excited 
s-like level of the KB-pseudopotential, and so there is a ghost state beneath the Cu 4s- 
state. By contrast, taking the s-component as the local potential results in a ghost-free 
KB-pseudopotential. 



3 Description of the package 

The package fhi98PP provides as main utilities the UNIX C-shell scripts psgen and pswatch. 
The script psgen and the related FORTRAN program fhipp solve the all-electron atom 
and generate the ionic pseudopotentials (see Sees. 2.1 - 2.3). The script pswatch and the 
related FORTRAN program psip solve the pseudo-atom for a given pseudopotential, and 
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carry out tests to monitor the pseudopotential's transferability in its semilocal and fully 
separable form (see Sees. 2.4 - 2.6). The shell scripts serve as command line interfaces for 
running the FORTRAN programs. For visual inspection of results, they group output into 
thematic graphic files to be viewed with the public domain, X- Windows based plotting 
tools XVGR or XMGR [69] . The scripts include command summaries (execute psgen -h 
or pswatch -h). The Tables 2-4 explain the required and the optional input of psgen 
and pswatch. The Tables 5-7 summarize the meaning of the input and output files, 
and display their dependencies on the FORTRAN routines. All input and output files 
reside in the chosen working directory. The name of the input files can be chosen freely. 
For the output files, the shell scripts compose filenames from a freely chosen mnemonic 
identifier and preset functional suffixes. For example, "al.cpi" shall denote an output file 
which tabulates an aluminum pseudopotential. 

The supplement "Test run protocol" at the end of this paper provides a prototype session 
protocol for generating pseudopotentials, for aluminum as an example. A database of input 
files is distributed along with the package (see Sec. 3.3). Below we outline the structure of 
the shell scripts and the related FORTRAN programs, and then discuss the various input 
variables. 

3.1 The tool psgen for solving the all-electron atom and generating pseudopotentials 

The script psgen reads from the commandline the identifier string name., and prinicpal 
input file <input data> containing the atomic configuration and run parameters as ex- 
plained below. 

psgen -o name <input data>, the standard mode, performs the all-electron calculation 
and constructs the pseudopotentials. Chief output is the run protocol name.dat and the 
pseudopotential data file name.cp\. 

psgen -o name -fc name. fc <input data>, the frozen-core mode, performs the all-electron 
calculation within the frozen-core approximation. There only the valence orbitals arc 
determined self-consistently while the core density is held fixed as given by the file 
name. fc (see below). This is meant to generate a certain pseudopotential component 
from an excited configuration with the unrelaxed core density of the ground state, and 
also to establish benchmark values, e.g. for excitation energies, in transferability tests, 
psgen -ao -o name <input data>, the atom-only mode, carries out just the all-electron 
calculation, as needed for instance to determine excitation energies. 

The option -e invokes an editor for the input file <input data> The command psgen -h 

lists all command options. 

The prinicpal input file {<input data>, e.g. called al.ini) specifies the electronic configu- 
ration of the atom, the partitioning of core and valence states, as well as parameters for 
handling exchange-correlation and constructing the pseudopotentials. Below we discuss 
the main variables, the Tables 2 and 3 give an overview. Table 4 describes a sample file. 
The variable z stands for the atomic number. Non-integer values to realize artificial atoms 
are allowed. 

The variables nc and nv specify the number of core states and the number of valence 
states, respectively. Note that the norm-conservation constraints can be met exactly for 
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one and only one valence state per angular momentum channel (see Sec. 2.2). 

The variable iexc specifies the exchange-correlation scheme used in solving the all-clcctron 
atom and in the unscreening of the pseudopotentials. Table 4 lists the implemented pa- 
rameterizations of the local-density approximation and the generalized-gradient approxi- 
mations (see Sec. 2.1). On the commandline iexc is set by psgen -xc iexc ... . 

The variable rnlc controls the handling of core-valence exchange-correlation in the un- 
screening of the pseudopotential and the construction of a partial core density (see Sec. 
2.3). Setting rnlc = corresponds to the ordinary linear treatment, a partial core density 
is not computed. Taking rnlc > implies the explicit treatment of core-valence XC with 
the help of a partial core density. The value of rnlc specifies the cutoff radius (in bohr) 
beyond which the partial core density is identical with the actual true core electron den- 
sity. The respective command is psgen -rc 1.3 e.g., setting rnlc — 1.3 bohr. The partial 
core density is appended to the file name.cpi. If the all-electron calculation is done in 
the frozen-core mode, the partial core density is constructed from the input core electron 
density. 

The variables n{i), l{i), and f{i) denote the i-th state's principal quantum number, angular 
momentum quantum number, and occupation number, respectively. The i = I, ...,nc core 
states need to be listed before the i — nc + 1, ...,nc + nv valence states. The principal 
quantum numbers should appear in ascending order, n{i) > n{i — 1). Unoccupied valence 
states which are not listed are assumed to be "unbound", the respective pseudopotential is 
then calculated by the Hamann procedure (see Sec. 2.2). Listing such states with f{i) = 
instructs the program to determine the respective state. Note that negatively charged 
ions, Zir=l™/(^) > ^) unstable within LDA or GGA, and may cause the calculation 
to fail. 

The variable Imax specifies the highest angular momentum channel up to which the 
program constructs pseudopotential components. The variable sjpp_def allows to select 
either the Hamann scheme s_pp_def = h or the Troullier-Martins scheme s_pp_def = t for 
constructing the pseudopotentials. Appropriate default cutoff radii are supplied by the 
program, psgen -T sjppJtype ... is the related script command. 

Any further input is optional and serves to override the default settings in the angular 
momentum channel It for the cutoff radius ret (in bohr), the reference energy ect (in 
eV), and the pseudopotential scheme sjppjtype. This feature is mainly for finetuning 
the cutoff radii when a ghost state of the Kleinman-Bylander form has to be eliminated 
(see Sec. 2.6). If the cutoff radii are enlarged in order to achieve softer pseudopotentials, 
the transferabihty of the resulting pseudopotential should be carefully verified. On the 
commandline, the options -rs, -rp, -rd, -rf may be used to set the cutoff radn of the s, p, d, f- 
components of the pseudopotentials. For instance psgen -rp 1.5 .. . sets the p-cutoff radius 
to 1.5 bohr. 

By the commandline option -r one instructs psgen to perform a non-relativistic all-electron 
calculation rather than a scalar-relativistic one (default). The corresponding variable tnrl 
is then set to true. 

By the commandline option -g all FORTRAN input and output files are saved, rather 
than being removed upon exiting psgen. 
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The script psgen first generates appropriate input files (fort. 20 and fort. 22) used by the 
FORTRAN routine fhipp which it invokes next. If the frozen-core mode is applied, the 
core density name.fc is copied to fort. 18 (see also below). After completion of fhipp the 
script writes a descriptive protocol to the file name.dat. It then generates the file name.cpi 
which contains the components of the ionic pseudotential, the pseudo wavef unctions, and, 
optionally, any partial core density on a logarithmic radial grid, adapted to the format 
used in the fhi96md electronic-structure and molecular- dynamics package [9]. If differ- 
ent reference configurations are used in constructing the pseudopotential components, a 
customized file name.cp\ can be created from the corresponding FORTRAN output files 
fort. 40, etc., described in Table 7. The output file name.aep contains the all-electron full 
potentials, and will be used by the utility pswatch to evaluate the all-electron logarithmic 
derivatives. The file name.fc lists the full core electron density, and may be used as in- 
put to psgen in a subsequent frozen-core calculation. Moreover the script produces various 
graphical files for visual inspection of, e.g., the ionic pseudopotential (file xv. name. pspotJ) 
or the all-electron wavefunctions (file xv. name. ae_wfct). 

In the pertinent FORTRAN program fhipp the main routine fhipp initially reads the run 
attributes and atomic configuration from the unformatted files fort. 20 and fort. 22. A log- 
arithmic radial mesh is set up, {r^ = <'^mc^h'"mm(^) I m = 1,2, ■■■,mjaax}, with default 
parameters as set in the FORTRAN header file default.h. All numerical integrations are 
done on this mesh. Any core electron density is read in from fort. 18, and must be tabu- 
lated on the current radial mesh. The subroutine sratom_n self-consistently solves for the 
all-clcctron atomic eigenstates and the full potential. The radial Schrodinger equations are 
integrated by a predictor-corrector scheme in the subroutine dftseq. The self-consistency 
cycle stops once all eigenvalues are converged. The program then writes a run protocol 
to fort.23, the all-electron wavefunctions to fort.38, the full core density to fort. 19, and 
the full all-electron potential to fort. 37. If the atom-only mode is apphed, the program 
exits at this point. Otherwise the main routine next determines the default parameters 
for constructing the pseudopotcntials, and checks for any optional input in fort. 22 to over- 
ride the default settings. The subroutine ncpp then determines the pseudo wavefunctions 
and the screened pseudopotcntials. The program proceeds with the unscreening of the 
pseudopotcntials. Optionally the subroutine dnicc constructs an appropriate partial core 
density, cither from the self-consistent full core density or from the user provided frozen- 
core density. Finally a protocol of the pseudopotential construction is appended to fort.23. 
The ionic pseudopotential and the pseudo wavefunction of each angular momentum chan- 
nel / are written to fort.4[/ = 0,1,...]. The components of the screened pseudopotential 
are written to fort.45 {I — 0), fort. 46 {I — 1), etc. . The partial core density is output to 
fort. 27. All quantities are tabulated on the dense logarithmic mesh that was employed in 
computing them, except for the all-electron wavefunctions which are output on a sparser 
mesh. 

3.2 The tool pswatch for solving the pseudo atom and checking pseudopotcntials 

The script pswatch is normally run after psgen. It reads from the command line the iden- 
tifying string name, and the name of the principal input file <input data> for the atomic 
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valence configuration. The latter has the same layout as when generating the pseudopoten- 
tials (sec Sec. 3.1). Furthermore pswatch assumes the presence of the files name.cp] (ionic 
pseudopotentials) and name.aep (all-electron potential). Both are routinely provided by a 
run of psgen as described above. The same radial grid must be used in these files. Also, the 
all-electron potential should be consistent with the present atomic configuration. Several 
run modes are implemented: 

pswatch -i name <input data>, the standard mode, performs first a self-consistent calcu- 
lation of the pseudo-atom, then computes the logarithmic derivatives for the all-electron 
potential and the screened pseudopotential, and carries out the ghost state analysis for 
the fully separable pseudopotential. 

pswatch -ao -i name <input data>, the atom-only mode, just solves the pseudo-atom, 
e.g., in order to determine excitation energies. 

pswatch -kb -i name <input data>, the semilocal-only mode, skips the ghost state anal- 
ysis. 

The option -e invokes an editor for the input file <input data>. The command pswatch 
-h gives a list of all options. 

Several parameters that control the evaluation of the logarithmic derivatives and the 
handling of the fully separable pseudopotential can be specified on the command line. 
The script pswatch collects them in an auxiliary input file (fort. 21). Below we introduce 
the main parameters, the Tables 2 and 3 list them in full. Sample input files are given in 
Table 4. 

The variable Imax specifies the angular momentum channel up to which pseudopotential 
components are present. For I < Imax the program uses the respective components to 
compute the logarithmic derivatives. For I > Imax it uses the local potential instead. 

The variable Hoc denotes the angular momentum channel whose pseudopotential com- 
ponent should be taken as the local potential, with Hoc < Imax. For instance the s- 
component is chosen by the pswatch -I ... command. 

The variable rlgd gives the diagnostic radius (in bohr) at which the logarithmic derivatives 
are evaluated. If rlgd = the program assigns as a default value the covalent radius for 
the atom at hand. Note that the diagnostic radius should be larger than the cutoff radii 
used in constructing the pseudopotential. The according command is, e.g., pswatch -rl 2.5 
. . . , setting rlgd — 2.5 bohr. 

The commandline option pswatch -cpo <outfile> ... causes the program to use for the anal- 
ysis of the Kleinman-Bylander form the self-consistently calculated pseudo wavefunctions 
rather than those provided by the file name. cpi (default). The calculated wavefunctions 
will be written to a file outfile analogous to name. cpi. This feature may be used, e.g., when 
just the ionic pseudopotential but not the respective pseudo wavefunctions are known. 
The corresponding variable tiwf true is then set true. 

The commandline option -r instructs pswatch to calculate the logarithmic derivatives for 
the all-electron potential non-relativistically rather than scalar-relativistically (default). 
The corresponding variable tnrl is then set true. 

By the commandline option -g all FORTRAN input and output files are saved, rather 
than being removed upon exiting psgen. 
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The script sets up appropriate input for the FORTRAN program psip which is invoked 
next. A self-evident protocol is printed to the terminal and also written to the file 
name.test. Upon completion, pswatch assembles the file xv. name. Ider for the visual in- 
spection of the logarithmic derivatives of the all-electron potential and pseudopotentials. 

In the FORTRAN program psIp the main routine psIp reads the run attributes and atomic 
valence configuration from the unformatted files fort. 21 and fort. 22, the latter having the 
same layout as for the fhipp program. Next it reads from fort. 31 the (logarithmic) ra- 
dial mesh, the pseudo wavefunctions to be used in the Kleinman-Bylander construction, 
the components of the ionic pseudopotential, and, if present, the partial core density. 
The program writes a run protocol to standard output. The subroutine psatom first self- 
consistently solves for the wavefunctions and the screening potential of the pseudo atom. 
The self-consistency cycle stops once all eigenvalues have converged. The program then 
writes the semilocal ionic pseudopotential and the computed pseudo wavefunctions to 
fort.4[/ = 0,1,...], the related screened pseudopotential to fort. 45, fort. 46, etc., and any 
partial core density to fort.27. If the atom-only mode is appUcd, the program stops at 
this point. Otherwise the subroutine ppcheck proceeds with the evaluation of logarithmic 
derivatives. This is done for a user-specified or a preset diagnostic radius. An adaptive 
energy mesh is employed, where the upper and lower bounds may be adjusted in the 
FORTRAN source file ppcheck.f. As the first step, the all-electron potential is read in 
from the file fort. 37. It must be given on the same radial mesh as the pseudopoten- 
tials. Integrating the radial Schrodinger equation for each angular momentum channel 
I, the all-electron logarithmic derivative is output to fort.5[Z = 0, 1, ...]. The correspond- 
ing logarithmic derivative of the self-consistently screened semilocal pseudopotential is 
tabulated in fort.6[/ = 0,1,...]. Next, ppcheck transforms the pseudopotential into the 
fully separable form (see Sec. 2.6). By default the inputted pseudo wavefunctions (file 
name.cpi) are used to set up the projector functions, optionally they are substituted by 
the pseudo wavefunctions obtained earlier in psatom. The subroutine derlkb then evaluates 
the logarithmic derivative for the respective nonlocal radial Schrodinger equation, which 
is output to fort.7[/ = 0, 1, ...]. In the next step, kibyii computes the bound state spectra 
for the semilocal and nonlocal pseudopotentials. Furthermore ppcheck examines the fully 
separable potentials for ghost states. Finally the subroutine kinkon evaluates the kinetic 
energy for each calculated pseudo wavefunction in momentum space (see Sec. 2.4). 

3.3 Set-up of the package 

The package fhi98PP is distributed as a tar-archive which can be unpacked by the UNIX 
command tar. The directory Dfhipp forms the root of the package's directory tree. 
In the directory Dpiay we supply sample input and output files for several exemplary pseu- 
dopotentials (aluminum, sodium, arsenic, selenium, and nitrogen), along with a detailed 
tutorial (as Postscript file). These examples should illustrate on a rather "hands-on" level 
the procedure for constructing pseudopotentials. 

The directory bin/Tools contains the shell scripts to generate (psgen) and validate 
(pswatch) the pseudopotentials. Also provided are some self-explanatory shell scripts 
(psdata_al, pstrans) to assist with the transferabihty tests discussed in Sec. 2.4. In the 
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directory bin/Elements we supply a database of inputfiles used in constructing the pseu- 
dopotentials. These have the generic names chemicaLsymboth .\n\ for the Hamann, and 
chemicaLsymbol:trr\.\n\ for the Trouhicr-Martins scheme. The directory bin/Xvgr contains 
the header files for the graphics output produced by the shell scripts. 

The directory src includes the sources and libraries of the program, as well as the Makefile 
for the compilation of the FORTRAN programs fhipp and psip that are invoked by the 
shell scripts. Both programs use linear algebra routines which need to be linked either 
from the ESSL-library or from a subset of the public domain LAPACK-routines. The 
latter is included as a library (libFREE) in the directory lib, and may be compiled either 
as specified in the Makefile, or, alternatively, simply along with the other routines. 

In order to run the package's tools psgen or pswatch one has to create the executables 
fhipp.x and psip.x. To this end one first needs to properly set in the Makefile the name of 
the FORTRAN compiler, the compiler and hnker fiags, and specify whether the ESSL or 
hbFREE hbrary should be invoked. Executing the UNIX command make all then does the 
compilation. Secondly, one needs to specify in the scripts psgen and pswatch the correct 
source paths for the FORTRAN executables and the graphics header files. We suggest 
to install the shell scripts psgen and pswatch as personalized commands. Note that for 
viewing the graphics output from psgen and pswatch one of the pubhc domain packages 
XMGR or XVGR must be available as well. 



4 Test run 



In the test run one calculates a Hamann-type pseudopotential for aluminum. To this 
end one executes psgen al -o al:h as command. This will copy the inputfile ahh.ini from 
the database bin/Elements into the present working directory, run the program fhipp and 
process its outputfiles into the various text and graphics files as described above. Secondly, 
one evaluates the logarithmic derivatives and performs the ghost state analysis. This is 
done by executing the command pswatch -i al:h ahh.ini, which will produce output on 
the terminal as well as to various text and graphics files described in Sec. 3. Samples of 
the principal output and a terminal session protocol are listed in the Section "Test run 
protocol" . A full set of the input and output files of the test run is given in the directory 
sample. Further examples can be found in the directory Dpiay. 
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Appendix 



This appendix lists total energy expressions computed by the utilities psgen and pswatch. 
All- electron atom - The total energy of the spherical atom for a given set of occupancies 
fni (corresponding to a ground- or excited state) is obtained from the radial wavefunctions 
Uni{r) and electron density p(r) = ^ T,ni fni\uni{r)\'^ as layed out in Sec. 2.1, 

^tot-AE[^] = T[p] + ^ ^H[^] ^ f--p(r) dh , (44) 

J r 

with the Hartree energy -E^[p] = ^ / V^[p] r]p{r) d^r (10), the exchange-correlation energy 
E^^ in LDA (12) or GGA (14), and the kinetic energy 

T[P\ = E fni^ni - I V[p; r]p{r) d'r , (45) 

nl 

stated in terms of the one-particle or orbital energy, J2ni fni^nu and the potential energy 
due to the effective potential V[p] r] (9). 

Frozen-core all-electron atom - When the atom is calculated in the frozen-core approxi- 
mation only the valence states are determined self-consistent ly. The core density is taken 
over from a previously solved atomic reference state, in the following indicated by zero 
sub- and superscripts. The effective potential (9) is obtained with the electron density 

p(r) = p^°'""(r) + p"'^' with (46) 



1 



,(0), 



1 



val 



PrW = 4^E/n/ </W ' P^-'\r) ^ -^^Y. fniKi{r)\' , (47) 



47rr^ 



nl nl 

and the kinetic energies of the core and valence electrons are given by 



cor6 

7^'^°"[Po] = E fni^ni - I VW, r]pr{r)d'r , (48) 

nl 



val „ 

T^^'M = E fniem - J V\p; ry^\r)d'r . (49) 

nl 

The atomic total energy in the frozen-core approximation then reads 

E'°'-^%o; p] = T--[po] + T-'[p] + E''^] + + / ~-p{r)d'r . (50) 

J r 

For the reference configuration, i.e. p = po, the unconstrained and frozen-core total ener- 
gies (44) and (50) take on the same value. 

Pseudo atom - For a given set of valence occupancies /; the pseudo valence states are 
determined by self-consistently solving the radial Schrodinger equations associated with 
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the semilocal pseudopotential, 



"_1 



(51) 



with the screening potential 



yHXC(^) ^ yXCj^ps ^ -core. ^] ^ yH[^ps. ^] 



(52) 



and pP**(r) = /; |M["^(r)| . The total energy of the pseudo atom is given by 



For the nonlinear core- valence exchange-correlation scheme, E-^'-" and V-^'^ depend on the 
particular partial core density pl°^^{r) (see Sec. 2.3). 
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Table 2 

Atomic configuration and pseudopotential parameters, given in the input file name.m of the 
scripts psgen and pswatch (file fort. 22 of the FORTRAN routines). 

variable type values 
required input, psgen and pswatch 

z real > atomic number (nuclear charge) 

nc integer > number of core states 

nv integer > number of valence states 

rnlc real = linearized core-valence exchange-correlation 

> nonlinear core-valence exchange-correlation, 

using a partial core inside the radius rnlc (in bohr) 

iexc integer exchanges-correlation'^ approximation 

1 LDA Wigner^ [49] 

2 LDA Hedin/Lundquist^ [48] 

3 LDA Perdew/Zunger^ '81 [47] 

4 GGA Perdew/Wang^'^ '91 [27] 

5 GGA Becke'^, Perdew^ '86 [30,31] 

6 GGA Perdew/Burke/Ernzerhof^'^ [28] 

7 LDA Perdew/Wang'^ '92 [29] + relativistic 

correction'^ of MacDonald/Vosko [59] 

8 LDA Perdew/Wang^ '92 [29] 

9 GGA Bccke^, Lec/Yang/Parr^ [30,32] 

10 GGA Perdew/Wangs '91, Lee/ Yang/Parr^ [27,32] 

n{i) integer principal quantum number of state i 

states should be in ascending order w.r.t. to n(i) 

l{i) integer angular momentum quantum number of state i 

f{i) real occupancy number of level i 

required input, psgen only 

Imax integer ... 4 maximum I — Zmax to generate pseudopotential for 

sjpp_def character h Hamann pseudopotential scheme 

t Troulher/ Martins pseudopotential scheme 

optional input, psgen only 

It integer ... 4 for It-th pseudopotential component: use alternate 

cutoff radius and/or reference energy and/or scheme 
from present line 
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variable 


type 


values 




ret 


real 


= 


core cutoff radius, use default 






> 


use this value (in bohr) instead 


et 


real 


= 


reference energy, use default: for an occupied state 








its eigenvalue, for an unoccupied state the highest 








eigenvalue of the occupied states 






7^ 


use this value (in eV) instead 


sjppjype 


character 




pscudopotential scheme, use sjpp-def 






h, t 


see required input 
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Table 3 

Upper section: Control parameters for the FORTRAN routine fhipp invoked by the script psgen, 
and passed by the input file fort. 20. Lower section: The same for the FORTRAN routine psip 
invoked by the script pswatch, and passed by the input file fort. 21. These files are created by 
and depend on the arguments given to the shell scripts. 

variable type values 

file fort. 20 of fhipp (optional input) 

tdopsp logical true compute pscudopotcntial (default) 

false atom-only mode: stop after solving all-electron atom 

tnrl logical true non-relativistic all-electron atom 

false dto. scalar-relativistic (default) 

file fort. 21 of psIp (required input) 

Hoc integer . . . ^max use lloc-th pseudopotential component as the local po- 

tential 

Ibeg integer . . . ^max lowest angular momentum channel to calculate logarith- 
mic derivatives for 

lend integer . . . ^max highest angular momentum channel to calculate logarith- 
mic derivatives for 

Imax integer Imax > angular momentum number up to which pseudopotential 

components are present for, apply the local potential in 
any higher channel 

rlgd real > diagnostic radius where the logarithmic derivatives are 

computed, ought be outside the core region 

tlgd logical true evaluate logarithmic derivatives (default) 

false dto. skipped 

tkb logical true analyze Kleinman-Bylander potentials (default) 

false dto. skipped 

tiwf logical true use input pseudo wavefunctions for the construction of 

Kleinman-Bylander potentials (default) 

false dto. but use calculated waves 

tnrl logical true assume non-relativistic all-electron calculation in the cal- 

culation of the logarithmic derivatives 

false dto. scalar-relativistic (default) 
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Table 4 

Sample input files to construct and validate an aluminum pseudopotential. 

Atomic configuration and pseudopotential parameters 

(file al.ini of shell scripts, file fort.22 of FORTRAN routines) 

13.00 3 2 8 0.0 z nc nv iexc rnlc 

1 2.00 n(i) l(i) f(i) 

2 2.00 

2 1 6.00 

3 2.00 
3 1 1.00 

2 h Imax s_pp_def 

1.25 0.00 h It ret et s_pp_type 

1 1.40 0.00 h 



control parameters of fhipp (file fort.20 created by psgen) 
.t. .t. tdopsp tnrl 



control parameters of psip (file fort. 21 created by pswatch) 

2 2 2 0.0 Hoc Ibeg lend Imax rid 

.t. .t. .t. .f. tlgd tkb tiwf tnrl 
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Table 5 



File dependencies for the psgen program. The upper part refers to input files , the lower one to 
output files. 





psgen files 


FORTRAN files and routines 


command op- 


input data 


fort.20 


fhipp 


tions 


t e ET v.fi.Tn.p inii 


fort.22 




fj.toTTiir ronfie''- 








uration 








core density 




fort. 18 


sratom_n 


fin TTnTiPTi-pnrp innnpi 








run protocol 


1 . 

namcdat 


fort. 23 


fhipp 








ncpp 


full core density 


name.fc 


fort. 18 


sratom_n 


full potential 


name.aep 


fort.37 


sratom_n 


ionic pseudopotentials, 


name.cpi 


fort.40 etc. 


ncpp 


pseudo wavefunctions, 




fort.27 


dnicc 


partial core density 








screened pseudopotentials 


XV. name. pspot_s 


fort.45 etc. 


ncpp 


ionic pseudopotentials 


XV. name. pspotJ 


fort.40 etc. 


ncpp 


pseudo wavefunctions vs. 


XV. name. ps_ae_wfct 


fort.39 


ncpp 


all-electron wavefunctions 








full core, partial core, 


xv.riame.density 


fort. 19 


sratom_n 


pseudo valence density 




fort.27 


ncpp 






fort.25 


dnIcc 


all-electron wavefunctions 


XV. name. ae_wfct 


fort.38 


sratom_n 
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Table 6 

File dependencies for the pswatch program. The upper part refers to input files, the lower one 
to output files. 



pswatch files FORTRAN files and routines 



command options, 


name.m 


fort.21 


psip 


atomic configuration 




fort.22 




ionic pseudopotentials 


name.cp'\ 


fort.31 


psip 


partial core density 


name.cpi 


fort.31 


psIp 


(if present) 








all-electron potential 


name.aep 


fort. 37 


ppcheck 


run protocol 


terminal, 


std.out 


psip 




name.test 




ppcheck 








kibyii 


ionic pseudopotentials, 


as specified 


fort.40 etc. 


psip 


calculated pseudo waves. 








partial core density 








pseudo wavefunctions 




fort.38 


psatom 


logarithmic derivatives, 


XV. name. Ider 


fort. 50 etc. 


ppcheck 


all-electron, semilocal. 




fort. 60 etc. 


derlkb 


fully separable case 




fort. 70 etc. 




ionic pseudopotentials 


XV. name. pspotJ 


fort.40 etc. 


psip 


screened pseudopotentials 


XV. name. pspotJ 


fort.45 etc. 


psip 
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Table 7 

Format of the pseudopotential file name.cpi and the related FORTRAN output files fort. 40, etc., 
and fort. 27. 



line 


col. 




C U J. ip 1 J.( Hi 


file name.cpi by script psgen : 




1 
1 


1 
1 


yion 


number of valence electrons 




2 




number of pseudopotential components 


2. ..11 






unused 


12 


1 


^max 


number of radial mesh points 




2 


C-mesh 


mesh increment Tm+i/Tm 


12... mmax + 11 


1 


m 


radial mesh index 




2 




radial coordinate (in bohr) 




3 




for / = 0: radial pseudo wavefunction. 






normalized as /q°° \v^^{r)\'^dr = 1 




4 


Vrirm) 


for Z = 0: ionic pseudopotential (in hartree) 


- for each Z = 1, 2, 


... ,1. 


nax : put a block like for Z = 


- if core- valence exchange-correlation applies, r"^'^ > 0: append fort. 27 


output fort.40, etc. by FORTRAN program 


fhipp : 


1 


1 


"#" 


marker 




2 


^max 


number of radial mesh points 




3 


'^mesh 


mesh increment rm+i/fm 




4 


/ 


current angular momentum channel 




5 


rf 


used cutoff radius 




6 


Z 


atomic number (nuclear charge) 


2 . . . Ulmax 1 


1 


m 


radial mesh index 




2 


I'm 


radial coordinate (in bohr) 




3 




radial pscudowavefunction 




4 




ionic pseudopotential 


output fort.27 by FORTRAN program fhipp 




1 . . . TTlmax 


1 




radial coordinate (in bohr) 




2 




partial core density, normalized as 

47r ff'^^^(r)r^ dr = iV'^™'', the circumstantial 

number of electrons in the partial core 




3 




dto. 1^* radial derivative 



dio. 2"'' rridial dcrivalivc 



39 



Test run protocol 



COMMANDLINE> psgen -v -o al:h al.ini 



13.00 3 2 8 0.0 



1 
2 
2 
3 
3 
2 

1 





1 


1 

h 

1. 

1. 



2. 
2. 
6. 
2. 
1. 

25 
40 



00 
00 
00 
00 
00 

0.00 
0.00 



z nc nv iexc rncl 
n 1 f 



h 
h 



Imax s_pp_def 

It ret et s_pp_type 



<QUIT> <RETURN> 



psgen - done: output 

al . ini 

al : h . aep 

al : h . cpi 

al : h . dat 

al:h.fc 

XV. al :h. ae_wf ct 
XV. al:h. density 
XV . al : h . ps_ae_wf ct 
xv.al:h.pspot_i 
xv.al:h.pspot_s 
XV . al : h . unscreen 

COMMANDLINE> more al:h.dat 

fhi pseudopotential tool fhipp - version rev 1998 



chemical symbol 
nuclear charge 
total charge 
number of core states 
number of valence states 
exchange-correlation model 
scalar-relativistic mode 
parameters radial mesh 



Al 

13.00 
.00 

3 
2 

8 LDA CA Perdew/Wang 1991 



493 



1 . 024700 



.480769E-03 



=== all-electron atom === 



n 



occupation eigenvalue (eV) 



40 



< 1 

< 2 

< 3 

< 4 

< 5 



1 
2 
2 
3 
3 





1 

1 



2.0000 
2.0000 
6.0000 
2.0000 
1.0000 



-1504.3103 
-107.5091 
-69.7240 
-7.8302 
-2.7840 



total energy 
kinetic energy- 
orbital energy 
coulomb energy 
hartree energy 
exchange-correlation energy 
xc potential energy 
number of iterations 
integrated density 

fhipp - all-electron atom done 

=== HAMANN mode === h 



(Hartree a.u.) 
-241.76605 
241.94185 
-134.51716 
-579.08935 
112.85767 
-17.47621 
-23.08499 
29 

13.00000 



convergence . OE+00 





1 


n 


radius : 


node 


peak 


default core 


X 





3 




.800 


2.023 


1.214 


X 


1 


3 




.800 


2.582 


1.549 


X 


2 


3 




.000 


.000 


1.549 



=== pseudo atom 
1 type rcore 



h 

1 h 

2 h 



1.2418974 
1.3692182 
1.5468790 



rmatch eigenvalue (eV) norm test slope test 

all-electron pseudo 1 = 1 = 

3.2163140 -7.8302025 -7.8302030 1.0000000 1.0000031 

4.4168841 -2.7839599 -2.7839599 1.0000000 1.0000024 

3.7233928 -2.7839599 -2.7839601 1.0000000 1.0000000 



linearized core-valence XC 

c-v equidensity radius 



total energy 
kinetic energy 
potential energy 
hartree energy 
xc energy 
integrated valence density 



1.43769 



(Hartree a.u.) 
-1.94588 
.62119 
-3.42557 
1 . 44497 
- . 58647 
3.00000 
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y range xvgr 



-3 11 



fhipp - done for input 






13.00 


3 


2 


8 .OOE+00 


z nc nv 





1 





2 


00 


n 1 f 


(§ 


2 





2 


00 




@ 


2 


1 


6 


00 




@ 


3 





2 


00 




@ 


3 


1 


1 


00 





2 

@ 1 



,25 
,40 



. OOE+00 
. OOE+00 



h 
h 



ImcLX s_pp_def 

It ret et s_pp_type 



COMMANDLINE> pswatch -v -i al:h al.ini 



2 2 2 0.0 
.t. .t. .t. 
13.00 3 2 



.f . 
8 



0.0 



1 

2 
2 
3 
3 
2 

1 







1 



1 

h 
1. 



2.00 
2.00 
6.00 



Hoc Ibeg lend Imax rlgd 
tlgd tkb tiwf tnrl 
z nc nv iexc rncl 
n 1 f 



2. 
1. 

25 



1.40 



00 
00 

0.00 

0.00 



h 
h 



Imax; s_pp_def 

It ret et s_pp_type 



<QUIT> <RETURN> 

fhi pseudopotential tool pslp - version rev 1998 



chemical symbol 
nuclear charge 
number of valence electrons 
number of valence states 
exchange-correlation model 
parameters radial mesh 
input pseudopotentials up to 



Al 

13.00 
3.00 
2 

8 LDA CA Perdew/Wang 1991 
493 1.024700 .480769E-03 
d 



=== pseudo atom (Hartree a.u.) === 



< n 1 occupation 

< 1 1 2.0000 

< 2 2 1 1.0000 



eigenvalue (eV) 
-7.8302 
-2 . 7840 



potential energy 
-1.21008 
-1.00541 
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total energy 


-1, 


. 94588 


kinetic energy 




.62119 


ionic pseudopotential energy 


-3, 


.42557 


hartree energy 


1, 


.44497 


xc energy 


-, 


. 58647 


local potential energy 


-3, 


.89662 


xc potential energy 




.76337 


integrated valence density 


3, 


.00000 


number of iterations 




16 


y range xvgr 




-3 1 



convergence 
1 



. OE+00 



pslp- pseudoatom done - now testing 

assuming scalar-relativistic all-electron atom 

d component taken as local potential 

input wavefunctions used for kb potentials 

kb potentials: spectrum of bound states (eV) 



1 eO el e2 

semilocal -7.8302 -.2102 .0000 

nonlocal -7.8302 -.2108 .0000 

semilocal 1 -2.7838 .0000 .0000 

nonlocal 1 -2.7838 .0000 .0000 



analysis of kb potentials: s waves 

* no ghost (ekb > 0, elocO < eref < elocl) 



kb cosine .3783 

kb energy 38.3109 eV ekb 

local potential groundstate -23.2494 eV elocO 

dto. 1st excited state -1.6612 eV elocl 

reference energy -7.8302 eV eref 



analysis of kb potentials: p waves 

* no ghost (ekb > 0, elocO < eref < elocl) 



kb cosine .3180 

kb energy 18.3317 eV ekb 

local potential groundstate -6.8226 eV elocO 

dto. 1st excited state -.0197 eV elocl 

reference energy -2.7840 eV eref 
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logarithmic derivatives: at radius = 2.9893 

nonlocal potentials 

all-electron potential 

semilocal potentials 



kinetic energy convergence in momentum space 





1 


n 


bracket 


cutoff 


norm 


kinet . energy 










(eV) 


(Ry) 




(Hartree) 


ck 







1 


. OE+00 


1 


.981164 


. 147167E+00 


ck 







1 


.OE-01 


9 


.999384 


. 179742E+00 


ck 







1 


. OE-02 


21 


.999971 


. 182613E+00 


ck 







1 


. OE-03 


30 


.999998 


. 182930E+00 


cx 





1 








1 . 000000 


. 182961E+00 


ck 


1 




1 


. OE+00 


2 


.959957 


.221304E+00 


ck 


1 




1 


.OE-01 


3 


. 998345 


.252144E+00 


ck 


1 




1 


. OE-02 


9 


. 999935 


. 254902E+00 


ck 


1 




1 


.OE-03 


16 


.999996 


.255235E+00 


cx 


1 


2 








1.000000 


. 255268E+00 



pswatch - done: output 
al :h. test 
XV . al : h . Ider 

COMMANDLINE> 
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